R Notebook

The goal here is to use conventional alpha diversity metrics to see how Chao1 richness, shannon diversity and evenness change across samples and to compare those to the values seen using breakaway in the AlphaDiversity.Rmd file

Setup

Run AlphaDiversity in scratchnotebooks That file calculates richness in breakawy which I will combine here

#source(here::here("RScripts", "InitialProcessing_3.R"))
source(here::here("RLibraries", "ChesapeakePersonalLibrary.R"))

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ───────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ ggplot2 3.3.6      ✔ purrr   0.3.4 
✔ tibble  3.1.8      ✔ dplyr   1.0.10
✔ tidyr   1.2.1      ✔ stringr 1.4.1 
✔ readr   2.1.3      ✔ forcats 0.5.2 ── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

ksource(here::here("ActiveNotebooks", "BreakawayAlphaDiversity.Rmd"))

processing file: /home/jacob/Projects/ChesapeaakeMainstemAnalysis/ActiveNotebooks/BreakawayAlphaDiversity.Rmd

  |                                                                                                                              
  |                                                                                                                        |   0%
  |                                                                                                                              
  |......                                                                                                                  |   5%
  |                                                                                                                              
  |...........                                                                                                             |  10%
  |                                                                                                                              
  |.................                                                                                                       |  14%
  |                                                                                                                              
  |.......................                                                                                                 |  19%
  |                                                                                                                              
  |.............................                                                                                           |  24%
  |                                                                                                                              
  |..................................                                                                                      |  29%
  |                                                                                                                              
  |........................................                                                                                |  33%
  |                                                                                                                              
  |..............................................                                                                          |  38%
  |                                                                                                                              
  |...................................................                                                                     |  43%
  |                                                                                                                              
  |.........................................................                                                               |  48%
  |                                                                                                                              
  |...............................................................                                                         |  52%
  |                                                                                                                              
  |.....................................................................                                                   |  57%
  |                                                                                                                              
  |..........................................................................                                              |  62%
  |                                                                                                                              
  |................................................................................                                        |  67%
  |                                                                                                                              
  |......................................................................................                                  |  71%
  |                                                                                                                              
  |...........................................................................................                             |  76%
  |                                                                                                                              
  |.................................................................................................                       |  81%
  |                                                                                                                              
  |.......................................................................................................                 |  86%
  |                                                                                                                              
  |.............................................................................................................           |  90%
  |                                                                                                                              
  |..................................................................................................................      |  95%
  |                                                                                                                              
  |........................................................................................................................| 100%

output file: /tmp/RtmpIdhQyV/file353597baadc

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio

Attaching package: ‘flextable’

The following object is masked from ‘package:purrr’:

    compose


Attaching package: ‘ftExtra’

The following object is masked from ‘package:flextable’:

    separate_header

Warning: Assuming taxa are rows

library(vegan)

Loading required package: permute
Loading required package: lattice
This is vegan 2.6-3

library(cowplot)
library(flextable)
library(ftExtra)

This file is dedicated to conventional, non div-net/breakaway stats, since breakaway seems to choke on this data.

Reshape back into an ASV matrix, but this time correcting for total abundance

nonSpikes %>% head

raDf <- nonSpikes %>% pivot_wider(id_cols = ID, names_from = ASV, values_from = RA, values_fill = 0)
raMat <- raDf %>% column_to_rownames("ID")

raMat1 <- raMat %>% as.matrix()

countMat <-  nonSpikes %>%
  pivot_wider(id_cols = ID, names_from = ASV, values_from = reads, values_fill = 0) %>%
  column_to_rownames("ID") %>% as.matrix()

seqDep <- countMat %>% apply(1, sum)
min(seqDep)

[1] 852

sampleRichness <- rarefy(countMat, min(seqDep))

rarefy everything to the minimum depth (852)

countRare <- rrarefy(countMat, min(seqDep))

Gamma diversity

specpool(countRare)

Doesn’t finish

#specpool(countMat)

Calculate diversity indeces

All richness estimates

richnessRare <- estimateR(countRare)

Shannon diversity

shan <- diversity(countRare)
shan

 3-1-B-0-2  3-1-B-1-2  3-1-B-180   3-1-B-20    3-1-B-5  3-1-B-500   3-1-B-53  3-1-S-0-2  3-1-S-1-2  3-1-S-180   3-1-S-20 
  4.358550   5.057168   4.645505   5.866334   5.195276   3.865637   5.506752   4.519812   4.961172   4.823132   5.314521 
   3-1-S-5  3-2-B-0-2  3-2-B-1-2  3-2-B-180   3-2-B-20    3-2-B-5  3-2-B-500   3-2-B-53  3-2-S-0-2  3-2-S-1-2  3-2-S-180 
  4.872777   4.372510   4.698438   4.575696   5.040413   5.437996   5.035825   4.843247   3.830293   4.950842   4.746431 
  3-2-S-20    3-2-S-5  3-2-S-500   3-2-S-53  3-3-B-0-2  3-3-B-1-2  3-3-B-180   3-3-B-20    3-3-B-5  3-3-B-500   3-3-B-53 
  4.889092   5.101738   4.806509   4.395303   4.241859   4.955895   3.457605   5.716900   5.220893   5.415232   5.436138 
 3-3-S-180   3-3-S-20  3-3-S-500   3-3-S-53  4-3-B-0-2  4-3-B-1-2  4-3-B-180   4-3-B-20    4-3-B-5  4-3-B-500   4-3-B-53 
  5.025812   4.927400   4.933217   4.319009   4.323957   4.946950   4.417338   4.456536   4.658738   4.345363   4.207244 
 4-3-O-1-2  4-3-O-180    4-3-O-5  4-3-O-500   4-3-O-53  4-3-S-0-2  4-3-S-180   4-3-S-20  4-3-S-500   4-3-S-53  5-1-S-1-2 
  5.033048   4.605692   5.147523   3.933170   4.501961   2.805419   4.558078   4.656982   4.669918   4.517247   4.362809 
 5-1-S-180   5-1-S-20    5-1-S-5  5-1-S-500   5-1-S-53  5-5-B-0-2  5-5-B-180    5-5-B-5  5-5-B-500   5-5-B-53  5-5-S-180 
  4.545527   4.163814   3.991503   4.054552   4.042100   4.632035   5.108552   5.484409   5.154378   4.903917   4.280018 
   5-5-S-5  5-5-S-500   5-5-S-53 C_5P1B_0P2 C_5P1B_180 C_5P1B_1P2  C_5P1B_20 C_5P1B_500  C_5P1B_53 
  4.934133   4.921967   4.306067   4.286760   4.825743   4.766759   5.454785   4.732364   5.196367 

Evenness

pielouJ <- shan/richnessRare["S.chao1",]
pielouJ

  3-1-B-0-2   3-1-B-1-2   3-1-B-180    3-1-B-20     3-1-B-5   3-1-B-500    3-1-B-53   3-1-S-0-2   3-1-S-1-2   3-1-S-180 
0.012650735 0.008180042 0.010628825 0.003710882 0.006572407 0.042372070 0.008906998 0.010251663 0.006769361 0.008460575 
   3-1-S-20     3-1-S-5   3-2-B-0-2   3-2-B-1-2   3-2-B-180    3-2-B-20     3-2-B-5   3-2-B-500    3-2-B-53   3-2-S-0-2 
0.007758092 0.008185653 0.013831141 0.010000060 0.012731121 0.006334166 0.006947808 0.007370699 0.010310991 0.017040735 
  3-2-S-1-2   3-2-S-180    3-2-S-20     3-2-S-5   3-2-S-500    3-2-S-53   3-3-B-0-2   3-3-B-1-2   3-3-B-180    3-3-B-20 
0.008120976 0.010779503 0.007686427 0.009669012 0.008608205 0.010571985 0.016670593 0.009355428 0.062019818 0.004974462 
    3-3-B-5   3-3-B-500    3-3-B-53   3-3-S-180    3-3-S-20   3-3-S-500    3-3-S-53   4-3-B-0-2   4-3-B-1-2   4-3-B-180 
0.007297796 0.006684979 0.005663804 0.007174607 0.009110878 0.011223098 0.013101420 0.009129840 0.007506505 0.009866906 
   4-3-B-20     4-3-B-5   4-3-B-500    4-3-B-53   4-3-O-1-2   4-3-O-180     4-3-O-5   4-3-O-500    4-3-O-53   4-3-S-0-2 
0.009604603 0.009595357 0.009468610 0.009368606 0.008100300 0.012180624 0.008613659 0.015215359 0.009862989 0.124685269 
  4-3-S-180    4-3-S-20   4-3-S-500    4-3-S-53   5-1-S-1-2   5-1-S-180    5-1-S-20     5-1-S-5   5-1-S-500    5-1-S-53 
0.010737086 0.007742281 0.012267734 0.011180515 0.010146067 0.011378040 0.015622114 0.016866247 0.011078011 0.014276923 
  5-5-B-0-2   5-5-B-180     5-5-B-5   5-5-B-500    5-5-B-53   5-5-S-180     5-5-S-5   5-5-S-500    5-5-S-53  C_5P1B_0P2 
0.008126376 0.008491069 0.006870717 0.009492710 0.007819712 0.013804485 0.011750911 0.009099696 0.016369217 0.006173593 
 C_5P1B_180  C_5P1B_1P2   C_5P1B_20  C_5P1B_500   C_5P1B_53 
0.007335071 0.007715700 0.003990014 0.007887274 0.004981453 

Combine diversity data

diversityData <- sampleData %>% 
  left_join(richnessRare %>% t() %>% as.data.frame() %>% rownames_to_column("ID"), by = "ID") %>%
  left_join(shan %>% enframe(name = "ID", value = "shannonH"), by = "ID") %>%
  left_join(pielouJ %>% enframe(name = "ID", value = "pielouJ"), by = "ID") %>%
  arrange(Size_Class)

Generate plots of diversity estimates

Parameters for all plots

plotSpecs <- list(
  facet_wrap(~Depth, ncol = 1) ,
  theme_bw(base_size = 16) ,
  geom_point(size = 4) ,
  geom_path(aes(color = as.factor(Station))) ,
  scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
  #scale_y_log10nice() ,
  scale_shape_manual(values = rep(21:25, 2)) ,
  scale_fill_viridis_d(option = "plasma") ,
  scale_color_viridis_d(option = "plasma") ,
  labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
  theme(legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = .5),
        axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)

Observed species counts, on rarefied data

plotObs <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = S.obs, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +ylab("Observed ASVs (Rarefied)")#+ scale_y_log10()
plotObs

Estemated Chao1 Richness

plotChao1 <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = S.chao1, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  geom_errorbar(aes(ymin = S.chao1 -2 * se.chao1, ymax = S.chao1 + 2* se.chao1), width = -.1) + 
  scale_y_log10() +
  ylab("Richness (Chao1)")
plotChao1

Shannon diversity

plotShan <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = shannonH, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  ylab("Diversity (Shannon H)") +
  lims(y = c(2.5, 6))
plotShan

Evenness

plotPielou <- diversityData %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +scale_y_log10() +ylab("Evenness (PielouJ)")
plotPielou

All plots together

plotAlpha <- plot_grid(plotObs, plotChao1, plotShan, plotPielou, nrow = 1, labels = LETTERS)
plotAlpha

ggsave(here::here("Figures", "ConventionalAlpha.png"), plotAlpha, width = 11, height = 4)

Do we see trends with lat and size?

Observed Species

Rarefied observed species numbers

obsMod <- lm(S.obs ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(obsMod)

Call:
lm(formula = S.obs ~ log(Size_Class) + I(log(Size_Class)^2) + 
    I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

Residuals:
     Min       1Q   Median       3Q      Max 
-223.262  -37.787   -0.657   47.789  195.661 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)          267625.988 133744.629   2.001 0.049325 *  
log(Size_Class)          32.892      8.103   4.059 0.000128 ***
I(log(Size_Class)^2)     -6.311      1.508  -4.186 8.24e-05 ***
lat                  -13937.556   6969.260  -2.000 0.049453 *  
I(lat^2)                181.530     90.726   2.001 0.049344 *  
depth                     4.420      3.209   1.377 0.172926    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 76.18 on 69 degrees of freedom
Multiple R-squared:  0.2636,    Adjusted R-squared:  0.2103 
F-statistic:  4.94 on 5 and 69 DF,  p-value: 0.0006353

Richness

Rarified chao1 estimates

chao1Mod <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(chao1Mod)

Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) + 
    I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

Residuals:
    Min      1Q  Median      3Q     Max 
-527.97 -117.95  -33.57  128.61  860.90 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)   
(Intercept)          842634.852 420822.680   2.002  0.04918 * 
log(Size_Class)          78.301     25.496   3.071  0.00305 **
I(log(Size_Class)^2)    -16.026      4.744  -3.378  0.00120 **
lat                  -43914.772  21928.525  -2.003  0.04915 * 
I(lat^2)                572.134    285.467   2.004  0.04898 * 
depth                    18.160     10.098   1.798  0.07650 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 239.7 on 69 degrees of freedom
Multiple R-squared:  0.1924,    Adjusted R-squared:  0.1339 
F-statistic: 3.289 on 5 and 69 DF,  p-value: 0.01011

As above but without latitude and depth

chao1ModSimple <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2), data = diversityData)
summary(chao1ModSimple)

Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2), 
    data = diversityData)

Residuals:
    Min      1Q  Median      3Q     Max 
-467.36 -149.31  -38.88  129.05  936.52 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           554.763     43.318  12.807  < 2e-16 ***
log(Size_Class)        78.966     25.726   3.070  0.00302 ** 
I(log(Size_Class)^2)  -16.380      4.785  -3.423  0.00103 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 242.2 on 72 degrees of freedom
Multiple R-squared:   0.14, Adjusted R-squared:  0.1161 
F-statistic: 5.859 on 2 and 72 DF,  p-value: 0.004389

Shannon Diversity

shanMod <- lm(shannonH ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

summary(shanMod)

Call:
lm(formula = shannonH ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityData)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.32065 -0.18411  0.01949  0.27550  0.85829 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           1.667e+03  7.858e+02   2.122   0.0375 *  
log(Size_Class)       2.024e-01  4.761e-02   4.250 6.57e-05 ***
I(log(Size_Class)^2) -3.673e-02  8.858e-03  -4.146 9.46e-05 ***
lat                  -8.661e+01  4.095e+01  -2.115   0.0380 *  
I(lat^2)              1.127e+00  5.331e-01   2.115   0.0381 *  
depth                 1.869e-02  1.886e-02   0.991   0.3250    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4476 on 69 degrees of freedom
Multiple R-squared:  0.3035,    Adjusted R-squared:  0.2531 
F-statistic: 6.014 on 5 and 69 DF,  p-value: 0.0001125

Evenness

pielouMod <- lm(pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(pielouMod)

Call:
lm(formula = pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityData)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.016294 -0.004747 -0.002583  0.000916  0.099924 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          -2.207e+01  2.632e+01  -0.839   0.4046  
log(Size_Class)      -3.952e-03  1.595e-03  -2.478   0.0157 *
I(log(Size_Class)^2)  6.757e-04  2.967e-04   2.277   0.0259 *
lat                   1.149e+00  1.372e+00   0.838   0.4049  
I(lat^2)             -1.494e-02  1.786e-02  -0.837   0.4055  
depth                -3.002e-04  6.316e-04  -0.475   0.6361  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01499 on 69 degrees of freedom
Multiple R-squared:  0.09657,   Adjusted R-squared:  0.0311 
F-statistic: 1.475 on 5 and 69 DF,  p-value: 0.2092

uomisto H (2010a). “A diversity of beta diver- sities: straightening up a concept gone awry. 1. Defining beta diversity as a function of alpha and gamma diversity.” Ecography, 33, 2–2

Prediction plots

Observed Species

predict(obsMod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12       13       14 
225.1538 225.1538 194.8219 194.8219 205.3509 156.7467 156.7467 184.4646 198.4497 300.2260 300.2260 269.8941 269.8941 280.4231 
      15       16       17       18       19       20       21       22       23       24       25       26       27       28 
231.8189 231.8189 259.5367 259.5367 331.0301 331.0301 300.6982 300.6982 311.2272 262.6230 262.6230 290.3408 304.3259 304.3259 
      29       30       31       32       33       34       35       36       37       38       39       40       41       42 
336.3394 336.3394 306.0076 306.0076 316.5366 316.5366 267.9323 267.9323 295.6502 295.6502 325.5520 295.2201 295.2201 305.7491 
      43       44       45       46       47       48       49       50       51       52       53       54       55       56 
305.7491 257.1449 257.1449 257.1449 284.8628 284.8628 298.8479 298.8479 295.0645 295.0645 264.7326 264.7326 275.2616 275.2616 
      57       58       59       60       61       62       63       64       65       66       67       68       69       70 
226.6574 226.6574 226.6574 254.3753 254.3753 268.3604 268.3604 255.1192 224.7873 224.7873 235.3163 235.3163 186.7121 186.7121 
      71       72       73       74       75 
186.7121 214.4299 214.4299 228.4151 228.4151 

$se.fit
 [1] 26.65044 26.65044 25.92384 25.92384 25.78549 27.43275 27.43275 28.97311 33.61118 18.69492 18.69492 18.08628 18.08628 16.98512
[15] 19.78897 19.78897 21.23350 21.23350 19.01263 19.01263 18.55660 18.55660 17.02002 19.83509 19.83509 20.93929 27.82376 27.82376
[29] 19.01048 19.01048 18.52392 18.52392 16.74918 16.74918 19.36508 19.36508 20.31747 20.31747 18.32932 17.69829 17.69829 15.79767
[43] 15.79767 18.21194 18.21194 18.21194 19.19052 19.19052 26.22172 26.22172 19.07012 19.07012 18.19426 18.19426 16.47139 16.47139
[57] 18.19575 18.19575 18.19575 19.25451 19.25451 25.78762 25.78762 24.15321 23.21205 23.21205 22.06503 22.06503 22.85093 22.85093
[71] 22.85093 23.83784 23.83784 28.86054 28.86054

$df
[1] 69

$residual.scale
[1] 76.17928

diversityData$pred_obs = predict(obsMod, se.fit = TRUE)$fit
diversityData$se_obs = predict(obsMod, se.fit = TRUE)$se.fit

plotSpecs2 <- list(
  facet_wrap(~Depth, ncol = 1) ,
  theme_bw(base_size = 16) ,
  #geom_point(size = 4) ,
  geom_path(aes(color = as.factor(Station))) ,
  scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
  #scale_y_log10nice() ,
  scale_shape_manual(values = rep(21:25, 2)) ,
  scale_fill_viridis_d(option = "plasma") ,
  scale_color_viridis_d(option = "plasma") ,
  labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
  theme(legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = .5),
        axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)

plotObs_pred <-  diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_obs, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_obs - 2 * se_obs, yend = pred_obs + 2 * se_obs, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5)  +
  plotSpecs2 + ylab("Predicted  ASVs") 
plotObs_pred

Richness

predict(chao1Mod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12       13       14 
461.6720 461.6720 360.3418 360.3418 441.7518 295.7073 295.7073 410.5479 376.5224 642.9478 642.9478 541.6177 541.6177 623.0276 
      15       16       17       18       19       20       21       22       23       24       25       26       27       28 
476.9831 476.9831 591.8238 591.8238 713.7117 713.7117 612.3816 612.3816 693.7915 547.7470 547.7470 662.5877 628.5621 628.5621 
      29       30       31       32       33       34       35       36       37       38       39       40       41       42 
719.9454 719.9454 618.6152 618.6152 700.0252 700.0252 553.9807 553.9807 668.8213 668.8213 687.4544 586.1243 586.1243 667.5342 
      43       44       45       46       47       48       49       50       51       52       53       54       55       56 
667.5342 521.4897 521.4897 521.4897 636.3303 636.3303 602.3048 602.3048 603.6376 603.6376 502.3074 502.3074 583.7174 583.7174 
      57       58       59       60       61       62       63       64       65       66       67       68       69       70 
437.6728 437.6728 437.6728 552.5135 552.5135 518.4879 518.4879 496.8539 395.5237 395.5237 476.9337 476.9337 330.8892 330.8892 
      71       72       73       74       75 
330.8892 445.7298 445.7298 411.7043 411.7043 

$se.fit
 [1]  83.85465  83.85465  81.56842  81.56842  81.13312  86.31617  86.31617  91.16284 105.75637  58.82290  58.82290  56.90783
[13]  56.90783  53.44308  62.26529  62.26529  66.81044  66.81044  59.82256  59.82256  58.38766  58.38766  53.55289  62.41039
[25]  62.41039  65.88474  87.54647  87.54647  59.81580  59.81580  58.28484  58.28484  52.70071  52.70071  60.93154  60.93154
[37]  63.92818  63.92818  57.67254  55.68702  55.68702  49.70679  49.70679  57.30320  57.30320  57.30320  60.38230  60.38230
[49]  82.50570  82.50570  60.00344  60.00344  57.24758  57.24758  51.82664  51.82664  57.25226  57.25226  57.25226  60.58363
[61]  60.58363  81.13983  81.13983  75.99721  73.03587  73.03587  69.42682  69.42682  71.89962  71.89962  71.89962  75.00489
[73]  75.00489  90.80864  90.80864

$df
[1] 69

$residual.scale
[1] 239.6954

diversityData$pred_chao1 = predict(chao1Mod, se.fit = TRUE)$fit
diversityData$se_chao1 = predict(chao1Mod, se.fit = TRUE)$se.fit

plotChao1_pred <-  diversityData %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_chao1, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_chao1 - 2 * se_chao1, yend = pred_chao1 + 2 * se_chao1, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5)  +
  plotSpecs2 + ylab("Predictd Richness (Chao1)") + scale_y_log10()
plotChao1_pred

Shannon Diversity

predict(shanMod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12       13       14 
4.500319 4.500319 4.326415 4.326415 4.297065 3.993768 3.993768 4.116157 4.383554 4.956811 4.956811 4.782907 4.782907 4.753557 
      15       16       17       18       19       20       21       22       23       24       25       26       27       28 
4.450260 4.450260 4.572649 4.572649 5.151681 5.151681 4.977777 4.977777 4.948427 4.645130 4.645130 4.767519 5.034915 5.034915 
      29       30       31       32       33       34       35       36       37       38       39       40       41       42 
5.197726 5.197726 5.023822 5.023822 4.994472 4.994472 4.691175 4.691175 4.813564 4.813564 5.145591 4.971687 4.971687 4.942337 
      43       44       45       46       47       48       49       50       51       52       53       54       55       56 
4.942337 4.639040 4.639040 4.639040 4.761429 4.761429 5.028825 5.028825 4.981512 4.981512 4.807608 4.807608 4.778258 4.778258 
      57       58       59       60       61       62       63       64       65       66       67       68       69       70 
4.474961 4.474961 4.474961 4.597350 4.597350 4.864747 4.864747 4.760193 4.586289 4.586289 4.556939 4.556939 4.253642 4.253642 
      71       72       73       74       75 
4.253642 4.376031 4.376031 4.643428 4.643428 

$se.fit
 [1] 0.15659117 0.15659117 0.15232183 0.15232183 0.15150895 0.16118783 0.16118783 0.17023857 0.19749069 0.10984658 0.10984658
[12] 0.10627035 0.10627035 0.09980025 0.11627494 0.11627494 0.12476260 0.12476260 0.11171336 0.11171336 0.10903381 0.10903381
[23] 0.10000530 0.11654591 0.11654591 0.12303394 0.16348531 0.16348531 0.11170074 0.11170074 0.10884180 0.10884180 0.09841393
[34] 0.09841393 0.11378429 0.11378429 0.11938023 0.11938023 0.10769839 0.10399060 0.10399060 0.09282305 0.09282305 0.10700868
[45] 0.10700868 0.10700868 0.11275861 0.11275861 0.15407212 0.15407212 0.11205113 0.11205113 0.10690481 0.10690481 0.09678168
[56] 0.09678168 0.10691355 0.10691355 0.10691355 0.11313459 0.11313459 0.15152148 0.15152148 0.14191809 0.13638804 0.13638804
[67] 0.12964847 0.12964847 0.13426621 0.13426621 0.13426621 0.14006503 0.14006503 0.16957713 0.16957713

$df
[1] 69

$residual.scale
[1] 0.44761

diversityData$pred_shanH = predict(shanMod, se.fit = TRUE)$fit
diversityData$se_shanH = predict(shanMod, se.fit = TRUE)$se.fit

plotShannonH_pred <- diversityData %>%

 filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_shanH, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_shanH - 2 * se_shanH, yend = pred_shanH + 2 * se_shanH, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"),  alpha = 0.5)  +
  plotSpecs2 + ylab("Predicted Diversity (Shannon H)") #+ scale_y_log10()
plotShannonH_pred

Evenness

predict(pielouMod, se.fit = TRUE)

$fit
          1           2           3           4           5           6           7           8           9          10 
0.019630086 0.019630086 0.021886750 0.021886750 0.021444103 0.024760944 0.024760944 0.022468058 0.019016295 0.010820847 
         11          12          13          14          15          16          17          18          19          20 
0.010820847 0.013077512 0.013077512 0.012634865 0.015951706 0.015951706 0.013658820 0.013658820 0.006908498 0.006908498 
         21          22          23          24          25          26          27          28          29          30 
0.009165163 0.009165163 0.008722516 0.012039357 0.012039357 0.009746471 0.006294707 0.006294707 0.005743555 0.005743555 
         31          32          33          34          35          36          37          38          39          40 
0.008000220 0.008000220 0.007557573 0.007557573 0.010874414 0.010874414 0.008581527 0.008581527 0.006479321 0.008735985 
         41          42          43          44          45          46          47          48          49          50 
0.008735985 0.008293338 0.008293338 0.011610179 0.011610179 0.011610179 0.009317293 0.009317293 0.005865530 0.005865530 
         51          52          53          54          55          56          57          58          59          60 
0.009217729 0.009217729 0.011474393 0.011474393 0.011031746 0.011031746 0.014348587 0.014348587 0.014348587 0.012055701 
         61          62          63          64          65          66          67          68          69          70 
0.012055701 0.008603938 0.008603938 0.013055314 0.015311978 0.015311978 0.014869331 0.014869331 0.018186172 0.018186172 
         71          72          73          74          75 
0.018186172 0.015893286 0.015893286 0.012441523 0.012441523 

$se.fit
 [1] 0.005244848 0.005244848 0.005101852 0.005101852 0.005074625 0.005398808 0.005398808 0.005701953 0.006614733 0.003679190
[11] 0.003679190 0.003559408 0.003559408 0.003342699 0.003894501 0.003894501 0.004178785 0.004178785 0.003741716 0.003741716
[21] 0.003651967 0.003651967 0.003349567 0.003903576 0.003903576 0.004120886 0.005475760 0.005475760 0.003741293 0.003741293
[31] 0.003645536 0.003645536 0.003296266 0.003296266 0.003811079 0.003811079 0.003998509 0.003998509 0.003607239 0.003483050
[41] 0.003483050 0.003109006 0.003109006 0.003584138 0.003584138 0.003584138 0.003776725 0.003776725 0.005160476 0.005160476
[51] 0.003753029 0.003753029 0.003580658 0.003580658 0.003241595 0.003241595 0.003580951 0.003580951 0.003580951 0.003789318
[61] 0.003789318 0.005075045 0.005075045 0.004753390 0.004568167 0.004568167 0.004342432 0.004342432 0.004497098 0.004497098
[71] 0.004497098 0.004691323 0.004691323 0.005679799 0.005679799

$df
[1] 69

$residual.scale
[1] 0.0149922

diversityData$pred_pielouJ = predict(pielouMod, se.fit = TRUE)$fit
diversityData$se_pielouJ = predict(pielouMod, se.fit = TRUE)$se.fit

plot_pielouJ_pred <- diversityData %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_pielouJ - 2 * se_pielouJ, yend = pred_pielouJ + 2 * se_pielouJ, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Evenness (Pielou J)") + scale_y_log10()
plot_pielouJ_pred

Combined prediction plot

plotPredictions <- plot_grid(plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred, nrow = 1, labels = LETTERS)

Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (geom_segment).

plotPredictions

ggsave(here::here("Figures", "ConventionalAlphaPredictions.png"), plotPredictions, width = 11, height = 4)

Even combindeder

plot_grid(plotObs, plotChao1, plotShan, plotPielou,
          plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred,
          nrow = 2, labels = LETTERS)

Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (geom_segment).

Combined summary table

alphaSummary <- tibble(
  metric = c("Observed ASVs", "Richness (Chao1)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
  model = list(obsMod, chao1Mod, shanMod, pielouMod)
)

alphaSummary <- alphaSummary %>%
  mutate(df = map(model, ~broom::tidy(summary(.))))

alphaSummary <- alphaSummary %>%
  select(-model) %>%
  unnest(df)

alphaSummary <- alphaSummary %>%
  rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
  mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
         Term = str_replace(Term, "\\^2", "\\^2\\^"),
         Term = str_replace(Term, "depth", "Depth"),
         Term = str_replace(Term, "lat", "Latitude"),
         Term = str_replace(Term, "_", " ")# BOOKMARK!!
         ) %>%
  mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
           reformat_sci()
         ) %>%
  mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
           reformat_sci()
  ) %>%
  mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
  mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
  rename(`Standard\nError` = `Standard Error`) %>%
  identity()

alphaSummary %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>%
  bold(i = ~ p< 0.05, j = "p") %>%
  colformat_md() %>%
  set_table_properties(layout = "autofit") %>%
  align(j = -c(1:2), align = "right")

Metric	Term	Estimate	Standard Error	T Value	p
Observed ASVs	Intercept	2.7 x 105	1.3 x 105	2.00	0.049
	log(Size Class)	3.3 x 101	8.1 x 100	4.06	< 0.001
	log(Size Class)2	-6.3 x 100	1.5 x 100	-4.19	< 0.001
	Latitude	-1.4 x 104	7.0 x 103	-2.00	0.049
	Latitude2	1.8 x 102	9.1 x 101	2.00	0.049
	Depth	4.4 x 100	3.2 x 100	1.38	0.173
Richness (Chao1)	Intercept	8.4 x 105	4.2 x 105	2.00	0.049
	log(Size Class)	7.8 x 101	2.5 x 101	3.07	0.003
	log(Size Class)2	-1.6 x 101	4.7 x 100	-3.38	0.001
	Latitude	-4.4 x 104	2.2 x 104	-2.00	0.049
	Latitude2	5.7 x 102	2.9 x 102	2.00	0.049
	Depth	1.8 x 101	1.0 x 101	1.80	0.077
Diversity (Shannon H)	Intercept	1.7 x 103	7.9 x 102	2.12	0.037
	log(Size Class)	2.0 x 10-1	4.8 x 10-2	4.25	< 0.001
	log(Size Class)2	-3.7 x 10-2	8.9 x 10-3	-4.15	< 0.001
	Latitude	-8.7 x 101	4.1 x 101	-2.11	0.038
	Latitude2	1.1 x 100	5.3 x 10-1	2.11	0.038
	Depth	1.9 x 10-2	1.9 x 10-2	0.99	0.325
Evenness (Pielou J)	Intercept	-2.2 x 101	2.6 x 101	-0.84	0.405
	log(Size Class)	-4.0 x 10-3	1.6 x 10-3	-2.48	0.016
	log(Size Class)2	6.8 x 10-4	3.0 x 10-4	2.28	0.026
	Latitude	1.1 x 100	1.4 x 100	0.84	0.405
	Latitude2	-1.5 x 10-2	1.8 x 10-2	-0.84	0.406
	Depth	-3.0 x 10-4	6.3 x 10-4	-0.48	0.636

Now considering breakaway values

richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.))))

Warning: `rename_()` was deprecated in dplyr 0.7.0.
Please use `rename()` instead.

diversityDataWB <- full_join(diversityData,
                             richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.)))),
                             by = c("ID" = "break_sample_names"), suffix = c("", "_break")) %>%
  mutate(pielouJ2 = shannonH/break_estimate) %>%
  identity()

diversityDataWB

pielouMod2 <- lm(pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityDataWB)
summary(pielouMod2)

Call:
lm(formula = pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityDataWB)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.013930 -0.005148 -0.002493  0.000898  0.105967 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          -1.957e+01  2.755e+01  -0.710   0.4799  
log(Size_Class)      -3.287e-03  1.669e-03  -1.969   0.0530 .
I(log(Size_Class)^2)  5.741e-04  3.106e-04   1.848   0.0688 .
lat                   1.018e+00  1.436e+00   0.709   0.4808  
I(lat^2)             -1.322e-02  1.869e-02  -0.707   0.4817  
depth                -2.380e-04  6.612e-04  -0.360   0.7200  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01569 on 69 degrees of freedom
Multiple R-squared:  0.06696,   Adjusted R-squared:  -0.000652 
F-statistic: 0.9904 on 5 and 69 DF,  p-value: 0.4301

Ok. So the narrative makes sense. Alpha diveristy is driven by variability in richness rather than evenness. Why would we see an effect in chao1 but not breakaway? Because chao1 is more driven by abundant stuff that makes the rarification threshold. My first hunch is that chao1 responds to evenness, but actually that shouldn’t have any effect since there is no evenness variability? Or maybe just that breakaway is more variable (because it detects fine level differences in rare species) and that doesn’t map as nicely with overall patterns.

plotBreak <- diversityDataWB %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  #scale_y_log10()+
  ylab("Richness (Breakaway)")
plotBreak

plotPielou2 <- diversityDataWB %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  #scale_y_log10()+
  ylab("Evenness (PielouJ)")
plotPielou2

Redo predictions for good measure

predict(pielouMod2, se.fit = TRUE)

$fit
            1             2             3             4             5             6             7             8             9 
 0.0116843027  0.0116843027  0.0135766875  0.0135766875  0.0133698701  0.0160079436  0.0160079436  0.0139748525  0.0099214190 
           10            11            12            13            14            15            16            17            18 
 0.0043268153  0.0043268153  0.0062192001  0.0062192001  0.0060123827  0.0086504562  0.0086504562  0.0066173651  0.0066173651 
           19            20            21            22            23            24            25            26            27 
 0.0011039449  0.0011039449  0.0029963297  0.0029963297  0.0027895123  0.0054275858  0.0054275858  0.0033944947 -0.0006589388 
           28            29            30            31            32            33            34            35            36 
-0.0006589388  0.0002124081  0.0002124081  0.0021047929  0.0021047929  0.0018979755  0.0018979755  0.0045360490  0.0045360490 
           37            38            39            40            41            42            43            44            45 
 0.0025029579  0.0025029579  0.0009065649  0.0027989497  0.0027989497  0.0025921323  0.0025921323  0.0052302058  0.0052302058 
           46            47            48            49            50            51            52            53            54 
 0.0052302058  0.0031971147  0.0031971147 -0.0008563188 -0.0008563188  0.0033197651  0.0033197651  0.0052121500  0.0052121500 
           55            56            57            58            59            60            61            62            63 
 0.0050053325  0.0050053325  0.0076434061  0.0076434061  0.0076434061  0.0056103150  0.0056103150  0.0015568815  0.0015568815 
           64            65            66            67            68            69            70            71            72 
 0.0066525974  0.0085449823  0.0085449823  0.0083381648  0.0083381648  0.0109762384  0.0109762384  0.0109762384  0.0089431473 
           73            74            75 
 0.0089431473  0.0048897138  0.0048897138 

$se.fit
 [1] 0.005490645 0.005490645 0.005340947 0.005340947 0.005312444 0.005651821 0.005651821 0.005969172 0.006924728 0.003851613
[11] 0.003851613 0.003726218 0.003726218 0.003499353 0.004077014 0.004077014 0.004374622 0.004374622 0.003917069 0.003917069
[21] 0.003823115 0.003823115 0.003506543 0.004086515 0.004086515 0.004314009 0.005732378 0.005732378 0.003916626 0.003916626
[31] 0.003816382 0.003816382 0.003450743 0.003450743 0.003989683 0.003989683 0.004185897 0.004185897 0.003776290 0.003646281
[41] 0.003646281 0.003254708 0.003254708 0.003752106 0.003752106 0.003752106 0.003953719 0.003953719 0.005402318 0.005402318
[51] 0.003928913 0.003928913 0.003748464 0.003748464 0.003393511 0.003393511 0.003748771 0.003748771 0.003748771 0.003966902
[61] 0.003966902 0.005312884 0.005312884 0.004976155 0.004782251 0.004782251 0.004545938 0.004545938 0.004707852 0.004707852
[71] 0.004707852 0.004911180 0.004911180 0.005945979 0.005945979

$df
[1] 69

$residual.scale
[1] 0.0156948

diversityDataWB$pred_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$fit
diversityDataWB$se_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$se.fit

plot_pielouJ2_pred <- diversityDataWB %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_pielouJ2 - 2 * se_pielouJ2, yend = pred_pielouJ2 + 2 * se_pielouJ2, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Evenness (Pielou J2)") #+ scale_y_log10()
plot_pielouJ2_pred

Breakaway richness subplots

plotBreakaway <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  geom_errorbar(aes(ymin = break_lower, ymax = break_upper), width = -.1) + 
  scale_y_log10() +
  ylab("Richness (Breakaway)")
plotBreakaway

#predict(breakLm, se.fit = TRUE)
# doesn't work because built with a different data frame

Why are these not smooth curves?!! What if I redo the model, this time with the same data frame

breakLm2 <- lm(break_estimate ~ log(Size_Class) + I(log(Size_Class) ^2) + lat +  I(lat^2) + depth ,data = diversityDataWB)
breakLm2 %>% summary()

Call:
lm(formula = break_estimate ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityDataWB)

Residuals:
    Min      1Q  Median      3Q     Max 
-2974.5 -1191.2  -151.6   599.9  6768.1 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          7124615.61 3339862.88   2.133   0.0365 *
log(Size_Class)          244.45     202.35   1.208   0.2312  
I(log(Size_Class)^2)     -75.16      37.65  -1.996   0.0498 *
lat                  -370568.38  174035.93  -2.129   0.0368 *
I(lat^2)                4817.28    2265.61   2.126   0.0371 *
depth                    151.10      80.15   1.885   0.0636 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1902 on 69 degrees of freedom
Multiple R-squared:  0.1414,    Adjusted R-squared:  0.0792 
F-statistic: 2.273 on 5 and 69 DF,  p-value: 0.0567

Note the non statistical significance overall

#predict(breakLm2, se.fit = TRUE)
diversityDataWB$pred_break = predict(breakLm2, se.fit = TRUE)$fit
diversityDataWB$se_break = predict(breakLm2, se.fit = TRUE)$se.fit

plot_break_pred <- diversityDataWB %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
#  filter(Station == 4.3) %>%
  ggplot(aes(x = Size_Class, y = pred_break, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_break - 2 * se_break, yend = pred_break + 2 * se_break, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Richness (Breakaway -- LM)") #+ scale_y_log10()
plot_break_pred

Rebuilding combined products

plotAlphaWB <- plot_grid(plotBreakaway, plotShan, plotPielou2, nrow = 1, labels = LETTERS)
plotAlphaWB

ggsave(here::here("Figures", "BreakawayAlpha.png"), plotAlpha, width = 11, height = 4)

Summary table I want both breakaway metrics here

bettaTable <- myBet$table %>% 
  as.data.frame() %>%
  rename(estimate = Estimates,
         `std.error` = `Standard Errors`,
         `p.value`=`p-values`
         ) %>%
  mutate(`statistic` = NA) %>%
  rownames_to_column(var = "term") %>%
  select(term, estimate, std.error, statistic, p.value) %>%
  as_tibble()
bettaTable

alphaSummary2 <- tibble(
  metric = c("Richness (Breakaway -- LM)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
  model = list(breakLm, shanMod, pielouMod2)
)
  
alphaSummary2 <- alphaSummary2 %>%
  mutate(df = map(model, ~broom::tidy(summary(.))))

## Add in willis variables

breakawaySummary <- tibble(
  metric = "Richness (Breakaway -- Betta)",
  model = NULL,
  df = list(bettaTable)
)

alphaSummary2 = bind_rows(breakawaySummary, alphaSummary2)

alphaSummary2 <- alphaSummary2 %>%
  select(-model) %>%
  unnest(df)

alphaSummary2 <- alphaSummary2 %>%
  rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
  mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
         Term = str_replace(Term, "\\^2", "\\^2\\^"),
         Term = str_replace(Term, "depth", "Depth"),
         Term = str_replace(Term, "lat", "Latitude"),
         Term = str_replace(Term, "_", " ")# BOOKMARK!!
         ) %>%
  mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
           reformat_sci()
         ) %>%
  mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
           reformat_sci()
  ) %>%
  mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
  mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
  rename(`Standard\nError` = `Standard Error`) %>%
  identity()



alphaSummary2

alphaTable2 <- alphaSummary2 %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>% bold(i = ~ p< 0.05, j = "p") %>% colformat_md() %>% set_table_properties(layout = "autofit") %>%
  align(j = -c(1:2), align = "right")
alphaTable2

Metric	Term	Estimate	Standard Error	T Value	p
Richness (Breakaway – Betta)	Intercept	7.1 x 106	2.4 x 102	NA	< 0.001
	log(Size Class)	1.2 x 102	6.1 x 101	NA	0.058
	log(Size Class)2	-5.0 x 101	1.2 x 101	NA	< 0.001
	Latitude	-3.7 x 105	6.1 x 100	NA	< 0.001
	Latitude2	4.8 x 103	1.6 x 10-1	NA	< 0.001
	Depth	1.5 x 102	1.0 x 101	NA	< 0.001
Richness (Breakaway – LM)	Intercept	7.1 x 106	3.3 x 106	2.13	0.036
	log(Size Class)	2.4 x 102	2.0 x 102	1.21	0.231
	log(Size Class)2	-7.5 x 101	3.8 x 101	-2.00	0.050
	Latitude	-3.7 x 105	1.7 x 105	-2.13	0.037
	Latitude2	4.8 x 103	2.3 x 103	2.13	0.037
	Depth	1.5 x 102	8.0 x 101	1.89	0.064
Diversity (Shannon H)	Intercept	1.7 x 103	7.9 x 102	2.12	0.037
	log(Size Class)	2.0 x 10-1	4.8 x 10-2	4.25	< 0.001
	log(Size Class)2	-3.7 x 10-2	8.9 x 10-3	-4.15	< 0.001
	Latitude	-8.7 x 101	4.1 x 101	-2.11	0.038
	Latitude2	1.1 x 100	5.3 x 10-1	2.11	0.038
	Depth	1.9 x 10-2	1.9 x 10-2	0.99	0.325
Evenness (Pielou J)	Intercept	-2.0 x 101	2.8 x 101	-0.71	0.480
	log(Size Class)	-3.3 x 10-3	1.7 x 10-3	-1.97	0.053
	log(Size Class)2	5.7 x 10-4	3.1 x 10-4	1.85	0.069
	Latitude	1.0 x 100	1.4 x 100	0.71	0.481
	Latitude2	-1.3 x 10-2	1.9 x 10-2	-0.71	0.482
	Depth	-2.4 x 10-4	6.6 x 10-4	-0.36	0.720

alphaTable2 %>% save_as_docx(path = here::here("Tables", "alphaTable2.docx"))

myBet$table

And finally predictions from richness, diversity evenness again.

plot_grid(plot_break_pred,plotShannonH_pred,plot_pielouJ2_pred, nrow = 1, labels = LETTERS)